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The temperature dependence of DNA flexibility is studied in the presence of stretching and 
unzipping forces. Two classes of models are considered. In one case the origin of elasticity is 
entropic due to the polymeric correlations, and in the other the double-stranded DNA is taken to 
have an intrinsic rigidity for bending. In both cases single strands are completely flexible. The 
change in the elastic constant for the flexible case due to thermally generated bubbles is obtained 
exactly. For the case of intrinsic rigidity, the elastic constant is found to be proportional to the 
square root of the bubble number fluctuation. 


I. INTRODUCTION 

To facilitate different fundamental biological processes, 
like replication, gene expression, assembly of functional 
nucleoprotein structures, and packaging of viral DNA, 
DNA has to go through a lot of twisting, stretching, and 
bending m- Generally different proteins induce these 
conformational changes in DNA, but not without fac¬ 
ing any resistance. This is because, when subjected to 
an external mechanical force, DNA responds elastically. 
Single-stranded DNA may be flexible and easy to bend, 
but double-stranded DNA (dsDNA) is known to be more 
rigid. However, the flexibility of short DNA fragments 
is important for different in vivo mechanisms, like those 
already mentioned, where loops or bends as short as 100 
base pairs in length are involved 0111 , and also in in vitro 
experiments, where fragments are used in open or hairpin 
geometries. It is therefore important to probe the elastic 
response of dsDNA not only in the thermodynamic limit 
of long length — dsDNA is long — but also for finite-size 
systems. 

Topological arguments, a la the Calugareneau theorem 
ttH, indicate the necessity of two major elastic constants 
of dsDNA, namely, the twist and the bending elastic con¬ 
stants. The former is tied to the helical nature of the 
double-helix and the latter is determined by both en¬ 
tropy and angular interactions between neighboring tan¬ 
gent vectors. These elastic moduli are characteristics of 
the bound phase and they disappear on DNAs melting 
into the denatured phase. It is quite analogous to the dis¬ 
appearance of the shear elastic modulus of a crystalline 
solid upon melting into the liquid phase. If dsDNA is 
treated as a free Gaussian polymer with noninteracting 
monomers, even then it exhibits an entropic elasticity 
originating from the correlations of a random walk. On 
the other hand, an intrinsic rigidity against bending at 
short scales (a semi flexible chain) produces a temper¬ 
ature dependent persistence length (~ 150 base pairs), 
within which a dsDNA acts more or less like a rigid rod. 
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Thus, it seems, a larger force is required to bend DNA 
of a length shorter than its persistence length than DNA 
of a longer length nmn]. Recent debates [T5H25] on 
the behavior of short segments brought into focus the 
importance of broken base pairs on its eventual or effec¬ 
tive rigidity. As the base-pair energy is comparable to 
the thermal energy at physiological temperatures (~2-3 
kcal/mol), bubbles form spontaneously or are produced 
by external forces (see, e.g., [2^ for an earlier study). 
Consequently, the issue of the elastic response of dsDNA 
cannot be studied in isolation as its intrinsic property 
but rather needs to be coupled to the inner degrees of 
freedom, namely, base pairings responsible for the bound 
state. 

Breaking the base pairings can separate dsDNA into 
two independent single strands and this melting hap¬ 
pens at a particular temperature. The thermal melting of 
DNA is by itself an interesting problem and important for 
different in vivo or in vitro processes. A notable example 
is the polymer chain reaction, which is used extensively 
in DNA amplification. Other than that it has been pro¬ 
posed [27U29] that at the dsDNA melting point the addi¬ 
tion of a third strand may support a three-stranded DNA 
bound state where no two pairs of strands are expected to 
be bound. This novel three-stranded DNA bound state 
is called Efimov-DNA and has been shown to support a 
renormalization-group limit cycle [301 EH- In the tem¬ 
perature region below the melting point there can be 
local melting at different positions, creating denatura- 
tion bubbles, which are nothing but single stranded loops 
preceded and succeeded by double-stranded segments. 
As single-stranded DNA is far more flexible than paired 
ones, these thermally generated bubbles can provide flex¬ 
ible hinges which can make dsDNA significantly flexible 
[EffiZl. Generally the average length of these bubbles 
increases as the critical point is approached and equals 
the system length at the melting temperature. The im¬ 
portance of bubbles for the melting transition is well un¬ 
derstood in the Poland-Scheraga framework [32] • The 
entropic contributions of the bubbles in different mod¬ 
els, originating from long range polymeric correlations of 
the individual strands, lead to different types of melting 
behavior, from weakly first order to infinite order EZl- 
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E]. The simpler coarse-gain level models show critical 
behavior [39l I42H44] . Close to melting, the bub¬ 

bles are therefore expected to contribute significantly to 
the flexibility, beyond just acting as hinges. Different 
from thermal melting is the unzipping transition, where 
the two strands of dsDNA are pulled apart at tempera¬ 
tures below the melting point. The unzipping transition 
is generally first order [151 Si], even in models with a 
first order melting transition EZIEHI. Since the unzip¬ 
ping force does not penetrate the bound state |42l [46] , 
the nature of the bubble distribution does not change in 
the presence of an unzipping force. As a result, the bub¬ 
bles are going to have their signatures on the flexibility of 
a DNA near the unzipping transition too. From a phase 
transition point of view, the bending elastic constant, de¬ 
spite its importance for DNA activity, is not a primary 
response function that characterizes a critical point. For 
example, one may compare it with the magnetic suscep¬ 
tibility of a ferro-para magnetic transition or the elastic 
modulus in a liquid-gas transition, whose divergence is 
associated with an exponent 7 . But as it is not the pri¬ 
mary response function, no such general results of critical 
phenomena are applicable here. Hence the necessity for 
a detailed study of the rigidity of melting DNA. In this 
paper we want to explore the elastic properties of DNA 
near its melting point. DNA melting is a genuine phase 
transition for which the DNA length has to satisfy the 
thermodynamic limit. But still, the existence of a tran¬ 
sition point is sufficient to affect a finite-size system even 
when it is away from the critical point. One of our aims 
is to obtain a few exact results on the elastic behavior 
for a class of models of DNA melting and unzipping. 

Different statistical mechanical models have been ap¬ 
plied with varying success to study the DNA elasticity 
problem. The classical semiflexible chain model with no 
denaturation bubbles has been employed by a number of 
investigators imiiH]. Segments made of single strands 
can be introduced in discretized semiflexible DNA, by 
considering models comprised of two-state internal coor¬ 
dinates and, also, by coupling these internal coordinates 
to the external rotational degrees of freedom of its tan¬ 
gent vectors [151 - I5T] . Bubbles appear naturally in our 
models without utilizing any other secondary variables. 

The modulus of interest comes from the response to a 
force applied at one end of each of the two strands, keep¬ 
ing the other end fixed at the origin. We use models of 
DNA where the strands are represented as polymers with 
native base pairing; namely, two monomers of the two 
strands interact only when they have the same contour 
length on the polymer. To probe the elastic behavior, we 
use a stretching force that would act on both strands in 
the same direction, while the phase can be changed by 
an unzipping force that acts on the strands in opposite 
directions. Here we quantitatively relate DNA flexibility 
to bubble-related quantities such as the bubble length, 
the bubble number etc. 

The organization of this paper is as follows. In Sec. [TT| 
we qualitatively describe the models, the flexible model 


and the rigid model, considered in this paper. Section 
hi] is devoted to the flexible model. In Sec. IIII A[ we in¬ 


troduce the corresponding recursion relations and define 
the observables necessary for analysis of both models. 
The elastic properties of the flexible model are explored 
in Sec. |HIB[ where we show a finite discontinuity in 
the elastic modulus at the melting point. We obtain the 
phase diagram in the presence of an unzipping force and a 
stretching force in Sec. HI D The transition is now first 
order and the elastic modulus shows a 5-function peak 
at the transition point. We introduce the rigid model in 
Sec. |IV[ After listing its governing recursion relations and 
defining the required observables specific to this model, 
we obtain its thermal melting (a continuous transition) 
point in Sec. |IV A| In Sec. |IVB[ we show that the cor¬ 
responding elastic modulus becomes anomalous around 
the melting point as it surpasses the unbound-state elas¬ 
tic modulus. The roles played by the bubbles are shown 
quantitatively in Sec. |IV C| Only the stretching force is 
considered for the rigid model. After a brief discussion 
of the relevance of our results in Sec. [^ we summarize 
and conclude in Sec. |VI| A few important supplementary 
materials are listed in the Appendices [A] and |B| 


II. QUALITATIVE DESCRIPTION 

To isolate the entropic and the intrinsic elastic behav¬ 
ior, two types of models are considered, viz., a flexible 
model, the standard model used for melting and unzip¬ 
ping [351 SOI I35H51] , and a rigid model, built from the 
flexible model. See Fig. [^ Both models show a zero- 
force melting point, generically denoted yc- The common 
features of these models are as follows : we consider each 
single strand as a directed polymer in a two-dimensional 
square lattice. We represent the base pairing by con¬ 
tact interactions between the monomers which can only 
occur at the same space and length coordinates. Cor¬ 
rect base-pair bonding is ensured by the directedness of 
the polymers. The chains are inextensible [58], of the 
same length, and attached to each other at the origin. 
We mimic DNA melting by the binding-unbinding phase 
transition in the system. The statistical weight of a con¬ 
tact interaction is the Boltzmann factor y = exp(,5e), —e 
being the energy per contact and j3 the inverse temper¬ 
ature with the Boltzmann constant set to fcs = 1. Such 
models in the past have been instrumental in studies of 
the melting and the unzipping transition of dsDNA and 
are known to show the relevant features of the higher di¬ 
mensional models [MESl nulls]- These models are also 
useful for studies of the dynamics of DNA. The flexible 
model [Fig. [^a)] has a hard-core repulsion that forbids 
the two chains to cross. The perfectly bound DNA re¬ 
mains as flexible as the single-stranded DNA so that at 
any nonzero temperature there is only the emergent en¬ 
tropic elasticity. In contrast, the rigid model [Fig. [^b)] 
has a bound state which has an intrinsic rigidity towards 
or against bending. In fact we consider the dsDNA to be 
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(a) 



FIG. 1. (color online). Schematic of dsDNA as two directed 
walks in the (1 + 1) dimension. The direction in which the 
forces act are indicated by arrows, (a) Flexible model : The 
polymers can not cross each other. The bound segments can 
bend left or right freely. This freedom can be restricted by 
introducing a statistical weight b. Here, b is associated with 
the left degree of freedom, (b) Rigid model : The polymers 
can cross here. The bound segments can not bend to the 
right and they are at least two bonds long, v is the weight 
associated with the bubble opening or closure. 

absolutely rigid in the bound state. The only way it can 
show any flexibility is via denaturation bubbles. Thus, 
the elastic response in this model is purely due to the flex¬ 
ible hinges made accessible by the bubbles. It is possible 
to allow some controlled semiflexibility, instead of abso¬ 
lute rigidity, via the introduction of a parameter 6, which 
penalizes the bound DNA taking an unfavorable turn. 
This model is discussed the Appendix This imposed 
rigidity is enough to give a melting transition even in the 
absence of any hard-core constraint. To incorporate this 
rigidity unambiguously, a constraint is required that a 
bound base pair can form only if its previous monomers 
are in the bound state. 

We apply an external space-independent mechanical 
stretching force independently at the end of each strand. 
If the two forces are in the same direction, the dsDNA is 
said to be under a stretching force gg. On the other hand, 
it will be under an unzipping force if the forces are in 
opposite directions. The average extension and the elas¬ 


tic modulus can be obtained from the free energy simply 
by taking derivatives with respect to the stretching force 
once and twice, respectively. 

We use the transfer matrix method through recursion 
relations to hnd the partition function of the system. For 
the analytical solution, the generating function for the 
grand partition function is used, from whose singular¬ 
ities the free energy can be determined [40j . For nu¬ 
merical calculations we iterate the recursion relations for 
finite lengths and find the exact partition function. The 
numerical calculations reflect the finite-size behavior of 
the concerning quantities. The effect of the unzipping 
force in the elasticity is also explored. The general case 
of two unequal forces can always be transformed into a 
case of unzipping and stretching forces. Since unzipping 
and stretching of dsDNA are independent of each other, 
we are able to generate a general phase diagram for three 
variables: the temperature, the stretching force, and the 
unzipping force. 


III. FLEXIBLE MODEL 

We use the model from Ref. [51] introduced to study 
the unzipping transition of dsDNA discussed in the last 
section. First, we solve the model analytically through 
the generating function technique, and then we study 
numerically the behavior of finite-length systems. 


A. Recursion relation and observables 

In the absence of any force the recursion relation fol¬ 
lowed by this system is given by |54j : 

Zn+lixi,X2)= ^ Zn{Xi+i,X2+j)[l-il-y)Sx^,x2], 

{i,j)=±l 

, ( 1 ) 

where Zn(xi, X 2 ) is the canonical partition function of the 
system of two polymers, each of length n, and the spatial 
positions of the nth monomers of polymer 1 and polymer 
2 are Xi and X 2 , respectively. For a given monomer num¬ 
ber if xi becomes equal to X 2 , then there is a contact. 
We set the initial condition as Zo(0, 0) = y such that 
two strands start from the origin. The non-crossing con¬ 
straint is implemented by not letting xi becoming greater 
than X 2 (xi < X 2 )- 

We apply a constant stretching force gs at the open end 
point of each strand. The partition function of n-length 
DNA in the presence of this stretching force is given by 

Z(5s) = X! ^n(xi,X 2 )e‘^‘'^, where X = Xi + X 2 (2) 

Xi ,X2 

and the sum is over all the allowed values of xi and X 2 - 
The elastic response of the system under a stretching 
force can be quantified through the average extension 
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X = 


df 1 dlnZjg,) 

dgs N dgs 


(x) and the elastic modulus (k). We define them in the 
following way 

dx 

and K = ——, (3) 

dgs 

where / = (3F = —\x\Z{gs) is the free energy of the 
system scaled by /?, and N is the length of the strands 
{N —>• oo). Using Eq. we can rewrite them as 

{X) 'Ex,,X2^Nixi,X2)eS‘‘^X 

X = — 


N 


^J2x,.X2^N{xi,X2)e9^^^ 




(4a) 


(4b) 


where (...) denotes the thermal average as indicated in 


Eq. (4a). An inspection of Eqs. (4a I and (4b) shows 


that X is related to the average vectorial position of the 
center of mass of the end points of the two strands of 
length n under a stretching force gs, and as expected, 
K is related to the fluctuation of x. If xi and X 2 are 
uncorrelated, then n is the sum of the individual elastic 
constants. This will be the case in the unbound phase. 
According to this definition for a given force the larger 
the value of k, the greater is the flexibility of the dsDNA. 
Two other important quantities are the average number 
of contacts between two strands {ric) and its fluctuation 
(Cc). Two extreme values of ric, 0 and 1, represent the 
unbound and the bound states, respectively. As ?/ is a 
temperature-like variable one can derive the specific heat 
of the system from Cc- We call it the specific heat for 
brevity. These are defined as 


y df duc 

and Cc = y^—. 
N oy oy 


(5) 


We follow these definitions in the rest of this paper. 

The recursion relations defined above can be solved ex¬ 
actly in the infinite-chain limit (i.e., the thermodynamic 
limit) with the help of generating functions. Based on 
the closed-form expressions, the physical quantities are 
obtained by taking appropriate derivatives as in Eq. ([^ . 
The results obtained in this way are called “analytical 
results.” These recursion relations can also be evaluated 
exactly, but numerically, by using the transfer matrix 
technique, for finite-length chains. The physical quanti¬ 
ties like the average extension and elastic modulus are 
then obtained by using Eqs. (|4a| and (|4b|. To evaluate 


Uc and Cc numerically we need to find out the first and 
the second numerical derivatives of the partition function 
with respect to y. We evaluate the recursion relations for 
the first and the second derivatives of the zero-force par¬ 
tition function by taking the first and second derivatives 
of Eq. ([^ with respect to y, respectively, and iterate 
them numerically to evaluate them. Then, using Eq. ([^ 
we calculate ric and Cc for zero applied force. For a con¬ 
stant nonzero applied force, we follow the same procedure 
but evaluate the partition function and its y derivatives 
by using Eq. (§ and taking derivatives of Eq. with 
respect to y. Such exact numerical results below are to 
called “numerical results.” 


B. Elastic Response Under a Stretching Force 



FIG. 2. (color online). Flexible model : Phase diagram of 
the system under a stretching force. The phase boundary 
separates the bound phase from the unbound phase. In the 
region 4/3 > i/ > 1 there exists a critical Qsc for every value of 
y . Above y — A the system is by default in the bound state. 
The solid blue line is the analytical curve, Eq. (11), and the 
red squares represent the nu meric ally obtained critical points 
[see discussion following Eq. (17)]. 


1. Generating function and the free energy 

By employing the generating function technique the 
recursion relation Eq. ([^ can be exactly solved. We 
define 

OO 

G{z,Xi,X 2) ='^z'^Zn{Xi,X2). ( 6 ) 

n—0 

By doing this we are going to the grand-canonical ensem¬ 
ble from the canonical ensemble. By multiplying both 
sides of Eq. 0 by and then summing over n we get 
two independent equations: one for nonzero unequal val¬ 
ues of Xi and X 2 , and another for Xi = X 2 = 0. These 
are given by 

-C{z,xi,X 2 )= C{z,xi + i,X 2 + j), (7a) 

— G(z,0,0) = --b C(z,i,j). (7b) 

yz z ^' 

(i,j)=±i 

The free energy per unit length of the DNA is determined 
by the singularity of C{z, xi,X 2 ) closest to the origin in 
the complex x-plane. 

Assuming a power-law form for G{z, xi,X 2 ) with re¬ 
spect to the relative position coordinate we make the 
ansatz 

G{z,xi,X 2 ) = (8) 
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where A and A are functions of z and independent of 
position coordinates. Equation ( ^ generalizes the ansatz 
in Ref. [M]. Using the ansatz, Eq. (Isl), in Eq. ([7a| and 


Eq. (7b I, two unknowns, A(gs, z) and \{gs, z), can easily 
be solved. Their forms are given in Appendix]^ The free 
energies of the different phases of the system are obtained 
from the singularities of G{gs,z). The singularity Z}, of 
A{gs,z) corresponds to the bound-state free energy and 
the branch point singularity z/ of X{gs,z) gives the free 
energy of the unbound state. They are calculated as 


Zbiy,9s) = 

zfiy^gs) = 


y- 1 


y sech{2gs) 
sech^(5s) 


/ sech^(2gs)~ ^ 

V y-1 


,(9a) 

(9b) 


The difference in the force term can be understood by 
looking at the low-energy excitations. In the case of the 
free chains a force gs flips a bond interchanging the en¬ 
ergies ±5s- This gives the sech^(5s) term. In the bound 
state, with coincident end points, a bound bond gets 
flipped under a force 2gs, yielding the sech^( 25 s) term. 
From here onwards we suppress y and ga as arguments for 
notational simplification and show them whenever neces¬ 
sary. The corresponding dimensionless free energies per 
unit length are given by 


/b = lnzb, (10a) 

//=lnz/. (10b) 


There are two parameters in this formulation, y and ga- 
The singularities move when these two parameters are 
changed. Consider a situation where the system is in the 
bound state with the free energy given by fb. Now, we 
can vary the parameters in such a way that the unbound- 
state singularity z/ crosses Zf, and becomes closest to the 
origin. In this situation the free energy of the system 
becomes //. The crossing of the singularities defines the 
transition point from bound to unbound by the force at 


( 11 ) 

The phase diagram in the y-gac plane is shown in Fig. 
The phase boundary has the following limiting forms: 

9sc ^ Vyc - y for y -b yc-, and (12a) 

^ for y-)-l + . (12b) 

iny 


2. Results and discussions 

a. Long-chain limit. When the two strands are in 
the unbound state, they come closer and form contacts 
in the influence of the stretching force. In this way an 
unbound state becomes a bound state above the critical 
stretching force. Figure|^shows how ric, calculated using 



FIG. 3. (color online). Flexible model : Variation of Uc with 
Qa for different values of Qu and y. A non-zero Uc indicates 
a bound state. For gu = 0, there is no unbound state for 
y> 4/3. The Qu = 0 curves are from Eqs. (|^, | |10a| |, and 
(lObl. The cyan line with triangles is for the melting point 
y ~ Vc and shows the quadratic dependence on ga- The black 
curve with filled diamonds represents y = 1.4 > y^. The red 
curve, for gu = 0, shows a continuous transition, while the 
blue curve, for nonzero gu, Eq. ( |21[ ), shows a discontinuity. 
The symbols on these analytically obtained curves, in the 
infinite-chain limit, are to make them distinct. 


Eq. becomes non-zero before saturating at 1 as the 
stretching force crosses a critical value for a y < yc- This 
shows that the transition is continuous. It is already 
known that at the point (4/3,0) in the phase diagram 
the system goes through a second-order phase transition. 
Beyond this critical point the system always remains in 
the bound state, thus excluding any other possibilities of 
phase transition. The only effect of the stretching force 
there is to influence the bubble statistics. The asymptotic 
behavior of Uc is given by 


iU9-9sc)Vy^c, ior 

for y = yc,ga 0, 


Uc 


where 


2 ' 

2 tis J 


no 

27 


^yviv-t) ’ 


[fiy-yc), 


no = ndy > yc,9s = 0 ) = 


for y > yc,ga -t 0, 
for y yc+,gs = 0, 

' - 2-b ^2/(2/- 1) 


2 ( 2 /- 1 ) 


(13) 


(14) 


That this is a second-order phase transition can be cor¬ 
roborated by examining the average extension of the cen¬ 
ter of mass due to the application of the stretching force 
and the elastic modulus. They are calculated using the 
definitions of Eq. Figure shows how the average 
extension changes continuously as the system crosses the 
critical stretching force. For a zero force x is 0, consis¬ 
tent with the Gaussian chain behavior, while the fully 
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FIG. 4. (color online). Flexible model : Plot of the average 
extension as a function of the stretching force with a fixed 
y = 1.20. These are obtained from Eqs. (^, ( |10a[ ), and (10b I. 
The average extension varies continuously around the critical 
point. The symbols on these analytically obtained curves are 
to make them distinct. 



FIG. 5. (color online). Flexible model: Elastic modulus curve 
for y = 1.20. The solid brown line is the analytically obtained 
curv e (see text). The dashed cyan line is for ki, = 4 sech^(2gs), 
Eq. (18aI, for the case with no bubbles. Other curves are the 
plots of K for different system lengths. Solid lines through the 
data points are guides for the eye. 


stretched state under a large force has x = 2. The slope 
discontinuity at the transition point gsc of Eq. 0 gives 
rise to a jump in the elastic constant as shown in Fig. 

To be noted here is that there is no pretransitional 
signature on either side of the transition. However, for 
a finite-size system the scenario is different. All other 
curves in Fig. [^except the analytical curve are for dif¬ 
ferent finite sizes of the system. In the unbound phase 
each strand has the equation of state x = tanh( 5 s) so 
that the total x = 2tanh((7s). This is the gs < gsc 


FIG. 6. (color online). Flexible model : A magnified version 
of Fig. around the same crossing point of all the curves. 
The chain length for each curve is given in the legend. Solid 
lines through the data points are guides for the eye. 



FIG. 7. (color online). Flexible model: The elastic modulus, 
Eq. (|15|), as a function of y for zero stretching force. 


branch. The corresponding entropic elastic constant is 
K = 2 sech^((/s). The completely bound state, in the ab¬ 
sence of any bubbles, should have a similar equation of 
state, with an elastic constant of purely entropic origin 
given by At = 4 sech^( 25 s). But the bubbles give an extra 
contribution. The exact form of the elastic constant can 
be determined for a few special cases. The y dependence 
of the zero force n is given by (see Fig. 


K(gs = 0 ) 



for y < yc, 
for y > yc- 


(15) 
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The elastic constant as a function of force at the melting 
point y = yc is 


, , 64w (w +1) , , 

K{y = yc) = — -with w = e^S\ (16) 

+ 14w + 1 )^' 


The behavior of n for y = yc and y > yc is shown in Fig. 



FIG. 8. (color online). Flexible model: The elastic modulus 
as a function of at j/ = 4/3 = j/c [Eq. ( |16| and y = 1.40 > yc 
[from Eqs. (|9a[) and (9b l]. 


b. Finite-length DNA. The contribution of the bub¬ 
bles becomes significant in finite-length DNA as shown 
in Fig. 1^ The finite-size effects become significant when 
the length is comparable to or shorter than the length 
of the bubble fluctuations. The elastic constant for a 
finite-length DNA is necessarily continuous, devoid of 
any singularity, but it should evolve into a discontinu¬ 
ous function as the length is increased. This indicates 
that shorter chains will show a larger deviation from the 
thermodynamic limit over a range of forces. A finite-size 
scaling form is 


= f((5s - 5sc)Af^^''), (17) 

so that K = f(0) at gg = gsc for all finite N. Therefore 
all the finite-size curves pass through a common point as 
shown in Fig. which is the critical point. By identify¬ 
ing the common points for other y values we can now de¬ 
termine the phase diagram numerically. This is shown in 
Fig-Il The consistency between this numerical method 
of finding the critical points with the analytical results 
helps us when the model under consideration is not solv¬ 
able analytically. All the points in this phase boundary 
including the thermal melting point (y = 4/3, = 0) 

are second-order critical points. The behavior of k shows 
that the system is most flexible when it is in the unbound 
state and under no external force, as it has the highest 
value of K. 


C. Role of the bubbles 


To highlight the importance of the bubbles we com¬ 
pare our results with the Y-model which is similar to 
the flexible model except that bubble formation is not 
allowed there [10]. The bound state of this model is the 
same as the completely bound state of the flexible model 
and it has a zero-force melting point (a first-order tran¬ 
sition) at yc =2. In the presence of gs the corresponding 
singularities and elastic constants are given by 


Zb 




1 

2y cosh 2gs ’ 
1 

2 cosh^ gs ’ 


Kfe = 4 sech^(2gs), and (18a) 
Kf = 2 sech^(gs), (18b) 


where Kb and Kf are the bound-state and the unbound- 
state elastic constants, respectively. We obtain the phase 
boundary by equating Zb with z/ as 


9sc — 


- cosh ^ 



(19) 


The phase boundary has similar asymptotics for y —>■ 
yci= 2) and ?/ —>■ 1 as in Eqs. (12aI and (12bI. In Fig. 
we compare Kb with the flexible model results. It shows 
that the flexibility of the bound state of the flexible model 
is mostly due to the bubbles. 


D. Elastic response in the presence of an unzipping 

force 


It is well known that this model undergoes an unzip¬ 
ping transition under the influence of an unzipping force 
in the absence of a stretching force. This unzipping tran¬ 
sition is known to be a first-order phase transition. The 
unzipped state consists of two completely separated in¬ 
dependent single strands. When DNA is in the double- 
stranded form the unzipping force tries to unzip it into 
two single strands. On the other hand, when DNA is 
in the unzipped state the stretching force tries to make 
them bound. Now, if we apply both the forces simulta¬ 
neously we expect a competition between the opposing 
effects. In this section we study this problem, again, ana¬ 
lytically for the infinite system and numerically for finite 
systems. We use the same definitions of quantities as in 
Eqs. (4a), (4b), and (§. 

Let us apply a spatially independent unzipping force 
gu at the rear end of the DNA, i.e., the forces act 
exactly in opposite directions. In the presence of a 
stretching force gs the generating function is given 
by G{gs,z) = A{gs, z)X{gs, where 

A{gs, z) and X{gs,z) are given by Eq. (Ala) and Eq. 
(Alb). The generating function in the presence of both 


forces is given by 


g{gs,gu,z)= ^ G(g«,z)eS“(—=). 

Xi ,X2 


( 20 ) 


















So, the bound-state singularity remains the same as 
Zb, Eq. (9a), consistent with the hypothesis of non¬ 
penetration of forces in the bound state [H], but the 
unbound-state singularity is now given by the solution 
of the equation = X{gs,z). This equation is easily 
obtained by performing the summation over xi and X 2 
in Eq. (20). Solving this equation for z we find that the 


unbound-state singularity is given by 




2 [cosh( 2 ( 7 s) -|- cosh(2g„)] ’ 


( 21 ) 


which corresponds to the partition function of the two 
chains under forces gs+gu and gs—gu, namely, 4cosh(5s-|- 
gu) cosh(gs —g„). For g^ = 0, the corresponding singular¬ 
ity matches Eq. (9b). The transition, as before, is given 


by the crossing of the singularities at 


guc = ^ cosh 


-1 


^ - cosh(2gs) 


( 22 ) 


This expression reduces to the known unzipping line m 
for gs = 0 and Eq. ([TT|) for g^ = 0. 



minima on the surface. The first-order surface ends on 
the ( 7 ii = 0 plane in a critical line that contains the usual 
melting point at ydgs = = 0). Except for the critical 

line all the other possible lines on the surface are first- 
order lines. To show that there is indeed a first-order 
transition we plot ric as a function of gs in Fig. with 
> 0 and y kept fixed. 



FIG. 10. (color online). Flexible model: x vs gs plot, with 
Qu = 0.4 and y = 1.20. Inset: Magnification of the critical 
region. The dotted magenta line is the analytical curve. See 
Sec. [m D 2[ All other curves are for finite system sizes shown 
in the legend. The discontinuity at Qsc indicates a first order 
transition. Solid lines through the data points are guides for 
the eye. 


2. Elastic constant 


FIG. 9. (color online). Flexible model: Three dimensional 
phase diagram of the system, Eq. (221, in the presence of 
both stretching and unzipping forces. All points on the sur¬ 
face except the curve for g^c = 0 represent first-order phase 
transition points. The unzipping line for = 0 is shown by 
the thick black line. 


1. Complete phase diagram with unzipping force 

After the introduction of an unzipping force we now 
have three control parameters. By changing any one of 
them while keeping the other two fixed, one can induce 
a phase transition in the system. The transition points 
are distributed on a surface in the y-gs-gu space given by 
Eq. (22). In Fig. [^we plot this function. The critical 
curve for guc = 0 in the surface is a second-order line. 
Around gs = 0, gudds) = 5«c(0) + ag'^ + ..., so that the 
unzipping line for gs = 0 lies along the locus of the local 


We plot X vs gs in Fig. [T0]and k vs gs in Fig. El keep¬ 
ing gu and y at fixed values. The dotted magenta curves 
are the analytical ones for an infinite system length, ob¬ 
tained by using Eqs. ( |^ , and ( |^ , while all the 


other curves are for finite system sizes which gradually 
match the analytical curve as N becomes larger, x shows 
a finite discontinuity at a critical gsc = 1.18. The ana¬ 
lytical curve for k has a d-function peak at gsc, which 
is not shown in Fig. E The uniform increase in the 
peak height with increasing system size in k at the crit¬ 
ical point is the signature of the delta peak. Below we 
list various useful limiting values of x and k. 


1. For gs-t 0,gu > gudy), 


x « 2 sech^ {gu)gs, 
K « 2 sech^ {gu) 2 


cosh(g^) - 2 
cosh'‘(gu) 


-9s- 


(23a) 

(23b) 
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1 

g, 


4. For 5s 5sc+, 5 = 1.2, 5 „ = 0.4, 

a; « 1.81211 + 0.66067(5s — ffsc), (26a) 
0.66067- 2.01486(5s-ffsc). (26b) 


For an infinite system the transition occurs suddenly at a 
single point. On the other hand, for a finite-size system 
the effect of the transition remains relevant for a domain 
of 5 s values containing 5sc beyond the scaling regime. 


IV. RIGID MODEL 


FIG. 11. (color online). Flexible model: k vs Qs plot with 
Qu = 0.4 and y = 1.20. Inset: Magnification of the critical 
region. The dotted magenta lines are the analytical cnrves. 
See Sec. |III D 2[ All other cnrves are for finite system sizes 
shown in the legend. The peak height increases proportionally 
with N, signaling a 5-function peak which is not shown in 
the analytical curve. Solid lines through the data points are 
guides for the eye. 


2. For 5s 0,5„ < 5„c(y), 


X 





8(3-25)v/^ 2 

y3/2 6's 


(24a) 

(24b) 


3. For 5 s 5 sc-, 2 /= 1.2,5„ = 0.4, 

a: « 1.57107-h0.73049(5s-5sc), (25a) 

K « 0.73049 - 1.03646(5s - 5sc). (25b) 


Here we customize the previous model to incorporate 
explicit weights for bubble formation. Doing that in the 
transfer matrix format is a bit involved. To identify a 
bubble we need to ensure that an unbound region is at¬ 
tached between two bound segments. A bound segment 
is defined as a DNA patch where every base pair is in 
the bound state and the minimum length it can have is 
2. We implement this by applying the constraint that 
a bound base pair can form only if another bound base- 
pair precedes it. So, for every step in the generation 
of the polymers we need to keep track of the previous 
step. We introduce a built-in rigidity to the dsDNA by 
instituting a bias against the bending towards the right 
of bound segments. For computational simplicity here 
we completely switch off the rightward option. By doing 
this we are introducing a bias in the propagation of DNA 
in favor of one direction. Other than the usual contact 
weight y we introduce another Boltzmann weight r; if a 
bound segment opens to form two single strands or two 
single strands recombine to form a bound segment. The 
recursion relation which obeys these rules is given by 


Zn{Xi,X2) 


'y[vZn-l{i,l) +VZn-lij,k) + Zn-l{j,l)], 
vz„-i{ij) +z„-i{i,k) + Zn-iij,k) -|-z„_i(j,/), 
vz„_i(5, k) + z„_i(j,fc) -I- Zn-l{i,l) + Z„-i{j,l), 
vyZn- 2 {i -I- 1,^ -I- 1) -I- z„-i{i, k) -I- z„-i(j, k) -|- z„_i(j, Z), 
vyzn-2{j -b l,fc -k 1) -k Zn-l{i,l) + Z„_i(t, fc) -k 2n-l(j,0: 

^Zn-l{i,k) + Zn-lii,l) +Zn-l{j,k) + Zn-l{j,l), 


if Xi = X2^n > 0 
if xi — X 2 = 2&n = 1 
if xi — X2 = —2&n = 1 
if xi — X 2 = 2&n > 2 
if Xi — X2 = —2&n > 2 
if \xi — X 2 \ > 2&n > 0 


(27) 


with i = xi — 1, j = xi + 1, k = X 2 — I and I = X 2 + I- 
The first two steps fix the initial configurations. To fix 
the configuration at the nth step we need to keep the 
information on not only the (n — l)th step but also the 
(n — 2)th step. We evaluate the recursion relation, Eq. 
(27), exactly for finite system sizes by iterating it nu¬ 
merically. Once Z(xi,X 2 ) is known, the force-dependent 


partition function can be obtained with the help of Eq. 
([^, and the corresponding elastic constant from Eq. 
Only the stretching force 5 s is considered here. Using the 
bubble weight v we can now count the average number 
of bubbles {nt) and calculate the average bubble length 
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(lb). They are formally given by the formulas 

df ^dnb (l-ric) 


where Cb describes the fluctuations in rif,. To evalu¬ 
ate them numerically, first we take the first and second 
derivatives of Eq. (27) with respect to v and iterate them 
numerically by setting u = 1 to obtain the first and the 
second derivatives of the zero-force partition function re¬ 
spectively. Once the derivatives of the partition function 
are known, Ub, Cb, and lb can be determined by using 
Eq. (28) for a zero applied force. To find out the cor¬ 


responding quantities for nonzero constant forces and to 
find out n, 

Section IIII Al 


and Cc, we follow the method described in 



1 1.2 1.4 1.6 1.8 2 


y 


A. Thermal Melting : = 0 


First, let us show that this model goes through a 
binding-unbinding transition as y is varied in the absence 
of any external forces. For the analysis in this section we 
set V = \. Here we use the exact numerical transfer ma¬ 
trix method for finite-size systems. 


In Fig. 12 we show how the average number of con¬ 
tacts varies with y. As the system size is increased one 
part of the curve gradually touches the y axis. And it is 
also evident that ric will saturate at ric = 1 for appropri¬ 
ately high y values. These indicate a binding-unbinding 
transition. To find out the order of the transition and 
the corresponding critical value of y, y^, we plot Cc vs y 
in Fig. Cc obeys a finite size scaling relation simi¬ 
lar to Eq. 0 which indicates a finite discontinuity. At 
y = 1.18 all the curves pass through a common point, 
implying yc = 1.18. The finite discontinuity at yc es¬ 
tablishes that this is a second-order phase transition. 
Note that we have not imposed the non-crossing con¬ 
straint here. The restriction imposed on the bound state 
is sufficient to induce a bound-unbound transition. In 
one spatial dimension the entropy of a system dominates 
over binding energy, which implies that there is no or¬ 
dering here. Imposition of special restrictions to limit 
the entropy may result in an energy dominated ordered 
state. The non-crossing constraint in the previous model 
did exactly that by decreasing the total number of con¬ 
figurations. The restriction imposed here on the degrees 
of freedom of the bound segment plays a similar role, de¬ 
creasing the total number of configurations of the DNA. 
We elaborate on this in Appendix [B| 


B. Elastic response ■. Qs ^ Q 

Let us now discuss the elastic properties of this sys¬ 
tem. As discussed earlier the inherent asymmetry in this 
model favors extension of DNA in one direction and op¬ 
poses it in the other direction. So under the influence of 


FIG. 12. (color online). Rigid model: ric increases from zero 
to non-zero values continuously at j/ > 1 before approaching 
the saturation value, 1. This indicates a binding-unbinding 
transition. Solid lines through the data points are guides for 
the eye. 



1 1.2 1.4 1.6 1.8 2 


y 

FIG. 13. (color online). Rigid model: Cc vs y plots for dif¬ 
ferent N. The peak height increases with increasing N but 
eventually saturates, creating a finite discontinuity. The dis¬ 
continuity occurs at ?/ = 1.18, where all the curves meet. Solid 
lines through the data points are guides for the eye. 


a spatially independent stretching force the DNA is more 
flexible in one direction compared to the other direction. 

As the system is no longer symmetric under o —gs 
we need to focus attention on the negative values of gs 
too. Figure [T4| shows that x varies continuously with gs, 
reaching ±2 for a large positive or negative gs- In Fig. 
[^we plot K vs gs keeping y fixed. Every curve shows a 
peak around gs = 0.02 which increases in height as N in¬ 
creases. The maximum of n, Kmax, goes to a finite value 
in the N ^ oo limit as shown in the inset in Fig. [T^ The 
inset in Fig. shows how Uc changes as we increase 
gs for y = 1.20. This indicates a continuous binding- 
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FIG. 14. (color online). Rigid model: Continuous stretching 
of the DNA to its full extent by both positive and negative 
forces. The y value is fixed at 1.20. Solid lines through the 
data points are guides for the eye. 



FIG. 15. (color online). Rigid model: k shows an increasing 
peak around Qs = 0.02 with increasing system length N and a 
fixed y = 1.20. Inset: Maximum values of k, Kmax, are plotted 
vs 1/A^. A linear fit with the fist four points (solid blue line) 
gives the estimate for N —>■ oo, Kmax = 4.95. This indicates 
that there is a finite discontinuity in k. The dashed black line 
is a plot of the function k = 2 sech^(gs), the unbound-state 
elastic constant. For Qs > 0.02, k matches the unbound-state 
modulus. Solid lines through the data points are guides for 
the eye. 

unbinding phase transition. Figure [T^ shows that Cc has 
a finite jump at gsc = 0.02 which can be identified as the 
common point in the peak region through which every 
finite-size curve passes. For < gsc, k shows anomalous 
behavior, as close to the transition point it can reach 
values which are much greater than the entropic elas¬ 
tic modulus of the unbound state given by 2 sech^( 5 s). 


shown in Fig. as a dashed black line. By collecting 



FIG. 16. (color online). Rigid model: Cc vs Qs curves for 
finite lengths as indicated. The y value is fixed at 1.20. Inset: 
Uc decreases continuously from a finite value to 0, indicating 
a continuous phase transition. Peak heights in Cc vs Qs curves 
saturate, indicating a finite discontinuity. Around Qs = 0.02 
all curves meet at the critical point. Solid lines through the 
data points are guides for the eye. 



FIG. 17. (color online) Rigid model: Numerical phase di¬ 
agram for the binding-unbinding transition. The solid line 
through the data points is guide for the eye. 

similar common points for different y values we draw the 
numerical phase diagram of the system, which is shown 
in Fig. [it] Noticeable there is that the stretching force 
actually unbinds the bound state for gs > gsc- The rea¬ 
son for this is the following. Due to the bias the bound- 
state formation on the positive x axis is unfavorable and 
the DNA prefers to go in the negative x direction. As 
gs is increased it wins over the bias eventually and pulls 
the DNA towards the positive x direction. Because the 
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bound state is forbidden in that direction the bound DNA 
unzips as a result. 


C. Role of the bubbles 

The flexibility of the bound state comes solely from 
the bubbles, as the bound segments are absolutely rigid 
in this model. Figure shows that rib becomes very 



-2 - 1.5 -1 - 0.5 0 0.5 1 1.5 2 

Ss 


FIG. 18. (color online). Rigid model: rib vs ga and lb vs Qs 
plot for y — 1.20. As the critical point is approached, rib 
becomes very small but 4 becomes as large as N. Solid lines 
through the data points are guides for the eye. 



FIG. 19. (color online). Rigid model: rib vs ga and Cb vs ga 
plot for y = 1.20. Around the critical point rib is very small 
but its fluctuation Cb is very large. Solid lines through the 
data points are guides for the eye. 

small while lb increases to almost equal to N as we ap¬ 
proach the transition point. Here lb « {N or 0) means the 
DNA is in the unzipped state. The peaks in rib curves in¬ 
dicate a large number of bubbles but at the same time of 



FIG. 20. (color online). Rigid model: ^y/rib vs y curves for 
different system lengths collapse to a single master curve in 
the bound region, flu = ffs = 0. flc = 1.18. Solid lines through 
the data points are guides for the eye. 



FIG. 21. (color online). Rigid model: vs y curves for 

different system lengths collapse to a single master curve in 
the bound region, flu = fls = 0. y^ = 1.18. Solid lines through 
the data points are guides for the eye. 


very small average length. From these two observations 
we can now say that as we approach the transition point 
many small bubbles coalesce to form large bubbles with 
decreasing numbers, and eventually becomes 0 when 
the two strands get completely separated. The fluctu¬ 
ation in rib, Cb, also becomes large around the critical 
point as shown in Fig. 19 Earlier we have shown that 
K in this model behaves anomalously and the anomalous 
behavior occurs in the same region where lb and Cb are 
the largest. In our model, n is the fluctuation of exten¬ 
sion X by definition and it depends on rib. For example, 
near the transition point rib is very small and x is also 
very small, although lb is large. This is because in the 
(1-1-1) dimension x for a single strand is 0 due to its 
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Gaussian nature. But for a bound DNA, x oc as 
shown in Fig. 20 It is then expected that k will be de¬ 


termined by the fluctuations in Ub, Cf,. In Fig. we 
plot vs y for different system sizes in the absence 

of any external force. For the bound region {y > 1.18) 
the curves collapse into a single master curve which is 
almost y independent, inferring that 


7.7x/a. 


(29) 


We therefore conclude that Cb is the important factor in 
determining the elastic behavior of the system. 


(1) For gs = gu = 0, 

k{9s = 0 ) = 


(2) For y = y^, i.e. at the critical point, 
&Aw {w + 1) 


for y < yc, 



Kiy = yc) = 


, with w = . 


-I- liw + 1)^^^ 

In the presence of the two opposing forces, gg as the 
stretching force and as the unzipping force, the tran¬ 
sition surface in the y-gs-Qu space is given by Eq. (22). 


For the model with intrinsic rigidity, the main result 
we obtained is 


7.7x/a. 


V. DISCUSSION AND SUMMARY 


There are a few points which we feel need to be clarified 
in more detail, (a) As the free ends of the two strands 
are stretched more and more with increasing forces in 
the same direction, they are bound to come closer due to 
their equal lengths and the starting ends’ being attached 
to each other. This coming closer together increases the 
possibility of the strands’ forming a bound base pair by 
gaining energy, (b) The transition in the presence of the 
unzipping force is definitely an unzipping transition. This 
is true even if gu is very small, (c) k in the flexible model 
is sensitive to the changes in gu, gs, and y. Elasticity 
is entropic in nature, which emerges from collective be¬ 
havior. The rigid model, on the other hand, has its own 
intrinsic elasticity. This is reflected in the anomalous be¬ 
havior of K for this model, (d) The single-molecule DNA 
experimental setup in which stretching force is achieved 
by placing the DNA in a directional flow can be a testing 
platform of our models, (e) In the nanopore sequenc¬ 
ing technique, dsDNA is unzipped and a single strand is 
passed through a nanopore [SI]. Other than that, during 
bacterial conjugation or infection of a cell by a virus the 
DNA assumes similar geometry. Our study may be rele¬ 
vant in these cases, (f) Single molecule DNA unzipping 
experiments are normally done at room temperature. In 
Eq. (22) we provided a phase boundary which depends 
not only on the temperature and the unzipping force, but 
also on the stretching force. Its remains a challenge to 
generate this phase boundary experimentally with tem¬ 
perature as a variable in single-molecule experiments. 


We summarize the basic results on rigidity as defined 
byEq. ^ for the two models. 

For the entropic rigidity as seen in the flexible-chain 
model of DNA, we obtained the following exact results. 


VI. CONCLUSION 

We have studied the effect of melting of DNA on 
its elasticity using (I - 1 - I) dimensional models by em¬ 
ploying exact numerical and analytical methods. Un¬ 
der a stretching force DNA goes through a second-order 
binding-unbinding phase transition. The dependence of 
DNA flexibility on the stretching force, the unzipping 
force, and the temperature has also been discussed. In 
the presence of both forces the system goes through a 
first-order unzipping transition. The complete phase di¬ 
agram in the y-gs-gu space is obtained. The average bub¬ 
ble length and the average bubble number for our model 
for different parameter values are also studied. We have 
shown that the DNA flexibility is related to the bub¬ 
ble number fluctuations. For zero external forces, the 
extension of the DNA is temperature independent and 
varies with the square root of the bubble numbers pro¬ 
portionally, while the elastic modulus is also proportional 
to the square root of the bubble number fluctuation. 
Though the binding-unbinding transition is very sharp 
for an infinite-length system, the transition point can in¬ 
fluence the elastic behavior of DNA for a broad region of 
parameter values when the system length is finite. Con¬ 
sequently, the elastic response of short-length DNA, as 
used extensively in experiments, has to be widely differ¬ 
ent from that of long-chain DNA. Furthermore, though 
DNA is a very long molecule, it can melt locally, depend¬ 
ing on the environment it is in. Thus our study will help 
us to understand the importance of these locally melted 
regions of shorter lengths in determining the elastic prop¬ 
erties of the DNA as a whole. 


Appendix A: A{gs,z) AND X{gs,z) 


Forms of A{gs,z) and X{gs,z) needed for Zb and Zf 
(Eqs. (j^ (9b)) are given by 
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A{9s,z) = -, ^ , 

cosh(2y,) ^(cosh(2y,) ^ 

(Ala) 

if 11 


K 9 s,z) = W 1 cosh( 25 y - ^1 + cosh(2gy. 

(Alb) 


Appendix B: BIAS-INDUCED MELTING 



1 1.2 1.4 1.6 1.8 2 
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FIG. 22. (color online), vs y plot for b = 0.5. ric becomes 
finite at y = 1.142. Numerical data are also consistent with 
the analytical result. Solid lines through the data points are 
guides are the eye. 


In the rigid model we have completely switched off one 
degree of freedom for the bound segments. We can do the 
entropy limiting job in a more general way by introducing 
a control parameter b instead, such that for b = 0 we 
block a degree of freedom of the bound state completely, 
for 5 = 1 we get back the good old free Gaussian chain, 
and for the intermediate values of b we obtain a partially 
biased scenario. The recursion relation followed by the 
system is now given by 

Z„+i(xi,xi) = y[Z„{xi + 1,311 + 1) + Znixi + l,a;i - 1) 

+ Znixi - l,Xi + 1) + bZnixi - l,Xi - 1)], 

Zn+l{xi,X2) = [Zn{xi + 1, X 2 + 1) + Zn{Xi + 1, X 2 — 1) 

+ Z„(xi — 1,X2 + 1) + Zn{xi — 1,X2 — 1)] 
for Xi ^ X2, (Bl) 


where .Z„(xi,X 2 ) is the partition function for the system 
length n. We set the initial condition as Zo(0,0) = y 
such that two strands start from the origin. The micro¬ 
scopic parameter 6 is a Boltzmann weight which controls 
the possibility of the two polymers’ both going to the 
right-hand side while being in the bound state. All the 
notations and definitions of the observables which are 
used in this model are the same as for the flexible model. 
By z-transforming Eq. (Bl| using the initial condition 



FIG. 23. (color online). Cc vs y plot for b = 0.5 showing a 
finite discontinuity at y = 1.142. Numerical data also show 
very similar behavior. Solid lines through the data points are 
guides for the eye. 


we get two independent equations for xi ^ X 2 ^ 0 and 

Xi = X2 = 0 : 

-G(z, xi, X 2 ) = G{z, xi -I- 1, X 2 -b 1) -I- G{z, xi — 1, X 2 -b 1) 

+G{z, Xi -b 1, X 2 — 1) 

+G{z, xi — 1, X 2 — 1), (B2a) 

— G{z, 0,0) = —bG(z,l,l) -bG(z,—1,1) 
yz z 

-bG(z, 1, —1) -b bGi^z, —1, —1). (B2b) 


To solve these independent equations we make an 
ansatz for the generating function, G(z,xi,X 2 ) = 
Substituting this ansatz in Eq. 
(B2a) and Eq. (|B2b) and solving for A and A we get 


-bz+ i + Z + 

2z -b yi - 4z - 1 
2z 


(B3a) 

(B3b) 


The bound-state and unbound-state singularities are 
given by 


6 - 1 - y 6 -b 1 -b yi/\/ 6 -b 1 - 46 -b 4 

=- (h \\2 -^P4a) 

( 0 - 



(B4b) 
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From these singularities all the other relevant quantities 
can be derived. Figure shows how Uc varies with y. At 
a high enough y, ric saturates to its maximum value, 1. 
In Fig. I^we plot Cc vs y. In both Fig. [2^ and Fig. 
we also show the corresponding numerical data for dif¬ 
ferent finite system sizes. Cc obeys a finite-size scaling 
form Cc = g{[y — yc)N^^^) such that y = yc = 1.142, 
Ccijjc) = 5 ( 0 )- The numerical data are consistent with 
the analytical results. From these two figures we con¬ 
clude that the system is going through the usual second- 
order binding-unbinding phase transition. We obtain the 
critical point of this transition analytically by matching 
Zb with Zf. The critical value of y is now b dependent 
and varies with b as yc = 4/(3 -I- 5). For 6 = 1 we have 
the same recursion relation as that of a system with two 
Gaussian chains which can freely cross each other and in 
which there is no phase transition. In our case also we 
get 2/c = 1, which means that there is no phase transition 
at hnite temperature. Another interesting limit is when 


^ = 0: 2 /c = 4/3, the critical point for two Gaussian chains 
with noncrossing constraint, although the recursion rela¬ 
tions for these two cases are not the same. Moreover, the 
b dependence of yc in our model adds extra flexibility in 
that we can now tune the critical point by tuning b for a 
wide range of values. 

From the recursion relation it is clear that the parame¬ 
ter & just modulates one of the four possible contributions 
to the partition function of the nth generation from the 
(n— l)th generation, b may be associated with any of the 
four possible contributions. The case when b is attached 
to the z{xi — I, a ;2 -I- 1) term is of special interest because 
the 6 —>■ 0 limit is exactly the flexible noncrossing case 
now. This current case also is exactly solvable through 
the generating function technique and gives exactly the 
same 6-dependent critical melting point, yc = 4/(3 -I- b). 
So, we can now actually modulate the noncrossing con¬ 
straint through b and the corresponding critical point as 
well. 
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